# -*- coding: utf-8 -*-
'''compute fir lowpass-filter coefficients

usage: python firdesign.py [arguments]

arguments:
    -h                 print this and exit
    -v                 be more verbose
    -F fN              Nyquist frequency (0.5), sets sps
    -C fCUT            cutoff-frequency (0.25)
    -s sps             sample rate (1), sets fN
    -W window          window name (hanning)
    -N coefficients    number of coefficients (100)
    -P                 plot frequency and time response functions
    -I filename        write coefficients as C code to filename

acknowledgements:
    Python code for basis FIR filter design
    @author: Matti Pastell <matti.pastell@helsinki.fi>
    http://mpastell.com
'''

import traceback
from pylab import *
import scipy.signal as signal

#Plot frequency and phase response
def plotresp(b, a = 1, fN = 0.5):
    figure()
    w,h = signal.freqz(b,a)
    h_dB = 20 * log10 (abs(h))
    subplot(211)
    f = w * fN / 3.141592654
    plot(f, h_dB)
    # ylim(-150, 5)
    ylabel('Magnitude (db)')
    xlabel(r'Frequency')
    title(r'Frequency response')
    subplot(212)
    h_Phase = unwrap(arctan2(imag(h),real(h)))
    plot(f, h_Phase)
    ylabel('Phase (radians)')
    xlabel(r'Frequency')
    title(r'Phase response')
    subplots_adjust(hspace=0.5)

    figure()
    l = len(b)
    impulse = repeat(0.,l); impulse[0] =1.
    x = arange(0,l)
    response = signal.lfilter(b,a,impulse)
    subplot(211)
    stem(x, response)
    ylabel('Amplitude')
    xlabel(r'n (samples)')
    title(r'Impulse response')
    subplot(212)
    step = cumsum(response)
    stem(x, step)
    ylabel('Amplitude')
    xlabel(r'n (samples)')
    title(r'Step response')
    subplots_adjust(hspace=0.5)

    show()


class Usage(Exception):
    def __init__(self, msg):
        self.msg = msg


class Error(Exception):
    def __init__(self, msg):
        self.msg = msg


def main(argv=None):
    import getopt

    if argv is None:
        argv = sys.argv

    options = "hvF:C:N:s:W:PI:"
    verbose = 0
    fN = 0.5
    fCUT = 0.25
    sps = 1.0
    window = 'hanning'
    N = 100
    do_plot = 0
    inc_file = None

    try:
        try:
            opts, args = getopt.getopt(argv[1:], options, ["help"])
            for opt, val in opts:
                if opt == '-h':
                    print __doc__
                    return 0
                elif opt == "-v":
                    verbose += 1
                elif opt == "-F":
                    fN = float(val)
                    sps = 2.0 * fN
                elif opt == "-s":
                    sps = float(val)
                    fN = sps / 2.0
                elif opt == "-C":
                    fCUT = float(val)
                elif opt == "-N":
                    N = int(val)
                elif opt == "-W":
                    window = val
                elif opt == "-P":
                    do_plot = 1
                elif opt == "-I":
                    inc_file = val

        except getopt.error, msg:
            print >> sys.stderr, __doc__
            raise Usage(msg)


        #Lowpass FIR filter
        #sps = 18.78
        #fN = 0.5 * sps
        #fcut = 2.5
        #n = int(3.0 * 18.78)
        #windowname = "hanning"
        a = signal.firwin(N, cutoff = fCUT, nyq = fN, window = window,
                          pass_zero = 1)

        print "n = %d   fcut = %.3f   fN = %.3f  win = %s" \
              % (N, fCUT, fN, window)

        if do_plot:
            plotresp(a, fN = fN)

        #print a

        if inc_file is not None:

            inc = file(inc_file, "w")
            print >> inc, "//"
            print >> inc, \
                  "// firdesign.py output for n = %d  fcut = %.3f  fN = %.3f" \
                  % (N, fCUT, fN)
            print >> inc, "// window = %s" % window
            print >> inc, "//"
            print >> inc, "const double coeff[] = {"
            for z in a:
                print >> inc, "%12.5e," % z
            print >> inc, "};"
            print >> inc, "const unsigned int ncoeff = " \
                  + "sizeof(coeff) / sizeof(coeff[0]);"
            print >> inc, "// end of firdesign.inc"


    except Usage, err:
        print >>sys.stderr, err.msg
        print >>sys.stderr, "for help use --help"
        return 2

    except RuntimeError, err:
        print >>sys.stderr, err
        return 3

    except:
        traceback.print_exc()
        return 1

if __name__ == "__main__":
    sys.exit(main())
    
